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Towards Arbitrary Accuracy Inviscid Surface 

Boundary Conditions 

Rodger W. Dyson* 

NASA Glenn Research Center at Lewis Field, Cleveland, OH 4-4135 

Ray Hixont 

Institute for Computational Mechanics in Propulsion, University of Toledo, Toledo, OH 43606 

Inviscid nonlinear surface boundary conditions are currently limited to 3 rd order ac- 
curacy in time for non-moving surfaces and actually reduce to order in time when the 
surfaces move. For steady-state calculations it may be possible to achieve higher accu- 
racy in space, but high accuracy in time is required for efficient simulation of multiscale 
unsteady phenomena. A surprisingly simple technique is shown here that can be used 
to correct the normal pressure derivatives of the flow at a surface on a Cartesian grid 
so that arbitrarily high order time accuracy is achieved in idealized cases. This work 
demonstrates that nonlinear high order time accuracy at a solid surface is possible and 
desirable, but it also shows that the current practice of only correcting the pressure is 
inadequate. 


I. Introduction 

H IGH order accuracy is required for efficient 
aeroacoustic simulations. 1-3 Recently, it has 
been demonstrated 4 ’ 5 that there is an optimal range 
of design accuracies in which the best efficiency and 
solution accuracy is obtained. The optimal value de- 
pends upon the length of the simulation and the level 
of accuracy required. For example, nonlinear invis- 
cid aeroacoutics simulations with 6 th order accuracy 
in space and time is optimal for propagating a wave 
five wavelengths with absolute error less than 10 -6 . 
Unfortunately, this is not typically achieved in prac- 
tice because of the solid boundary conditions 6 which 
may have high order spatial accuracy, but rarely ex- 
ceed 2 nd order time accuracy. This severely limits the 
efficiency of time marching methods. 

This paper demonstrates why most surface bound- 
ary approaches are typically only 1 st order accurate in 
time and it presents a surprisingly simple procedure for 
correcting the normal pressure derivatives to achieve 
any order of time accuracy in special cases. The ba- 
sic idea follows the work of Tarn 7,8 in which walls 
are modeled by pressure corrections; and follows the 
work of Hixon 9 in which the pressure corrections are 
treated separately from flow variables to greatly sim- 
plify their nonlinear implementation; and finally, those 
contributions are combined with Goodrich’s boundary 
treatment 10 to produce arbitrary time accuracy at a 


•Aerospace Engineer, Acoustics Branch, Structures and 
Acoustics Division, Research and Technology Directorate, 
t Senior Research Associate, Member AIAA 
Copyright © 2002 by the American Institute of Aeronautics and 
Astronautics, Inc. No copyright is asserted in the United States 
under Title IT, U.S. Code. The U.S. Government has a royalty- 
free license to exercise all rights under the copyright claimed herein 
for Governmental Purposes. All other rights are reserved by the 
copyright owner. 


surface. 

The usefulness of these ideas are currently limited to 
cases in which the pressure corrections are much larger 
than the density and velocity corrections. Nonetheless, 
the procedures shown here demonstrate the feasibility 
and efficiency of increasing the time accuracy at sur- 
faces. 


II. Problem Description 

The following form of the two-dimensional nonlinear 
Euler equations are used: 


dp dp dp 

M +U di + % + P 

dp dp dp 

m + '‘di + v fy +lp 


du dv 
dx + dy 
du 


dx 


+ 


d{pu) + d(pu 2 + p) + d(puv) 


dt 

d(pv) 

dt 


+ 


dx 
d(puv) 
dx 


dy 

d(pv 2 + p) 
dy 


Qi 

Qi 

Qz 

<?4 


(1) 

( 2 ) 

(3) 

(4) 


At the solid boundary we also have the following 
infinite set of boundary conditions: 10 


d a (V-n) 

dt a 


d a Q 5 
dt a 


Va : (a > 0) 


(5) 


Ordinarily, the source terms (Qi, Q 2 , Q 3 , Qi, QO 
are zero. But, for testing purposes, these terms will 
be modified so that the following analytical solution is 
defined: 5 


p(x, y , t) = ai cos(k x Trx) cos(k y Try ) cos(fc t 7rf) + Ci (6) 
u(x, y, t ) = o 2 cos(k x nx) cos(k y iry) cos(k t nt) + c 2 (7) 
v(x, y, t) = a 3 cos^Trx) cos(k y ny) cos(fc t 7rf) + c 3 (8) 
p(x, y, t) — a 4 cos(k x nx) cos(k y iry) cos(k t n t) + c 4 (9) 
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Due to the large number of derivatives used in what 
follows, the following shorthand notation will be used: 


x f _ d a+b+k f{x,y,t ) 

a,b,k Q x aybfk 


( 10 ) 


where x, when used, will indicate: updated (U), new 
(AT), old ( 0 ), or correction (C). 


III. Maintaining High Accuracy at 
Boundary 

This work uses high accuracy Hermitian Modified 
Expansion Solution Approximation (MESA) methods 
because high time accuracy at a surface requires cor- 
recting the spatial cross derivatives of the flow vari- 
ables, many of which the MESA scheme explicitly in- 
cludes on the grid. Any method could be used as long 
as the effect of cross derivative corrections are con- 
sistently applied. In Fig. 1, the absolute error of the 
numerical wavenumber for the derivative of a harmonic 
function, f(x ) = exp** 1 , is plotted by points per wave- 
length (PPW) as done in Ref. u The extraordinary 
resolution of the two-point MESA schemes allows one 
to use fewer grid points and larger time steps, resulting 
in significant computational efficiency for problems re- 
quiring high physical accuracy as shown in Fig. I. 5 ’ 15 
However, maintaining high accuracy at surfaces, 12 par- 
ticularly when waves are being scattered, 13 is required 
to maintain this high efficiency. 

Unfortunately, most inviscid surface boundary con- 
ditions are at most 2 nd order accurate in time since 
only the following conditions: 


d a (V -n) 
dt a 


= OVa : (a = 0, 1) 


( 11 ) 


are used when in fact Eq. (5) is true for any a. 

As an example, the numerical time advancement of 
the pressure on a surface expressed as a Taylor series 
is: 


p{x,y,t + At) = (12) 

Co, 0,0 + ^ 0 , 0 , 1 ^ + ^0,0,2 + • • • 

And from the governing equations Eqs. (1, 2, 3, 4) 
we can convert the series’ time derivatives into space 
derivatives. For example: 


“Co, o,i — 

Co,o, oCf,o, o + CoVoC^.o + 7Co,o,o (CVVo + CoYo) 

(13) 



a) MESA scheme compares well to popular techniques 



3rd Order 

5th Order 

7th Order 

9th Order 

11th Order 

13th Order 


15th Order 

17th Order 

- 1 9th Order 

21 st Order 


b) Resolution of Two-Point MESA schemes 



5 10 15 20 

-Log 10 (Error Level Required) 
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c) Efficiency of Nonlinear Euler Solver 

Pig. 1 






and 


r^p — 

°0,0,2 “ 

£*1,0,0 (£*0,0,0' £*?, 0,0 + £*0,0,0 £*0,1,0 + £*0,0^0 £*i,o,o) 

£* 0 , 0,0 (£ f , 0,0 £* 1 , 0,0 + £* 0 , 0 , 0 £* 2 , 0,0 + ^ 1 , 0 , 0 ^ 0 , 1,0 + 


>-'1,0,0 ^ L '0,0,0*-'l,0,0 "f °0,0,0'-'0,l,0 f '-'0,0, 0^1, 0,1 
Co, 0,0 (C“,o,oCf o t o + Co,o,oC 2 ,o,o + Ci,o,oCo,i,( 
Co.o.oCf lj0 + 7Cf,o,o (C“ 0 ,o + c 0 ” a>0 ) + 

7£*0,0,0 (£*2,0,0 + £*l,l,o)) + 

Co, 1,0 (Co,o,oCi,o,o + Co,o,oCo,i,o + ^o.oq^o.i.o) 
nv mu nv + Co i0i0 Cf jlo + Co >li0 Co il>0 + 
( /~ 1 U I /~iv \ i 



These equations clearly show that for 2 nd order time 
accuracy we must have the correct first order spatial 
derivatives of the pressure and velocities at the sur- 
face. In particular, Eq. (5) with a = 1 provides a 
constraint on these 1 st order spatial derivatives by con- 
verting the time derivatives to space derivatives again 
using Eqs. (31, 31): 


~Vx 

~Vy 


S r-iu r^u I nv nu 

°0,0,0 U 1,0,0 + ^0,0,0°0,1,0 

Co,o,oCi,o,o + Co, 0 ,oCo,i,o 


+ C 0 ,o^oCf, 0i0 l 
+ Co,o^oCo,i,oj 



which if not true, implies the time accuracy is only 1 st 
order. 


Similarly, the boundary condition Eq. (5) with a = 
2 will provide a constraint on the first and second or- 
der spatial derivatives (including density), which if not 
true, will reduce the numerical method to below 3 rd 
order time accuracy. And in general, one must enforce 
Eq. (5) for all a up to the desired time accuracy. Note, 
however, in the special case of a nonmoving inviscid 
wall that the even alpha conditions will automatically 
be satisfied if the initial flow conditions and their spa- 
tial derivatives are correct at the surface. 


A. Time Advancing Boundary 

If the spatial derivatives are correctly specified to 
maintain the desired time accuracy, then after time 
advancing the flow variables we find that the boundary 
condition, Eq. (5), with a — 0 will still be satisfied as 
shown below (via its Taylor series expansion in time): 


n x p(x, y,t + At)u(x, y,t + At)+ 
n y p(x, y,t + At)v(x, y,t + At) = 

= n x p{ x, y, t)u(x, y,t) + n y p(x, y, t)v(x, y,t)+ 



= p \n x u{x, y, t) + n y v(x,y, t)) + 

N — ^ - * 


o 


E 


i 

fc! 


=o 


(to: 

\ d k ~ l p 

j dt k ~‘ l 

v 

d l u 

Ux W 

A 

O l V 

V 



/ 

=0 ' 


(A t) k 


= 0 

(16) 

At the new time step, the boundary condition, 
Eq. (5), for a > 0 will in general not be true. For 
example, we know that after time advancing the prim- 
itive variables that for a = 1: 


n x {p(x, y, t + At)u(x, y,t + At)) t + 
n y (p(x, y,t + A t)v(x, y, t + A t)) t = 

= -n x ( p x u 2 + 2 puu x + Px + PyUV + pUyV + pUVy ) 
-Tly ( PyV 2 + 2 pVVy +Py+ P X UV + pU X V + PUV X ) 
+ 0 

(17) 

Clearly it is necessary to correct the spatial deriva- 
tives at each time step to maintain time accuracy at 
the boundary. In this work however, only the spatial 
derivatives of the pressure are corrected because this 
approach was successful for the 1 st order corrections 
used by Tam 7 ’ 14 and Hixon. 9 This approach reduces 
the seemingly intractible complexity of very high or- 
der nonlinear boundary conditions to a few lines of 
code. The downside to this approach is entropy and 
vorticity waves may be excited as no mechanism for 
properly correcting the density and velocity derivatives 
is provided here. For example, the one-dimensional 
characteristic form (entropy and two acoustic waves) 
of the Euler equations is: 


^ + A,^ = 0 for i = 1,2,3 
d(v 2 ) = d(u) + 

d(v 3 ) = d{u) - ^ ( 18 ) 

Ai = u, A 2 =u + c, A z—u-c 
vi = s — const, for dx = udt 
v 2 = u + — const, for dx = (u + a)dt 

v 3 = u — — const, for dx = (u — a)dt 

where c is the speed of sound, s is the entropy, d(v{) 
is the i th characteristic variable differential which is 
integrable for isentropic flows, and A, is the speed of 
the characteristic wave . 

The time derivative of the Euler equations provides 
a new set of characteristics with the same velocities 
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but include higher order space derivatives: 


(&)+(> 


(-A, 


) + A. 


+ Ai&® =0,i = l,2,3 


% dx 


(-^W) = (Ai)* 


d(ui) = d (-u{p x + £§-)) 


dVj 

dx 


d(v 2 ) = d l-(u + c)(u x + ^)) 
d(v 3 ) = d (-(u - c)(u x - ^)J 
Ai — u, A 2 = u + c, A 3 — u — c 


(19) 

where c is the speed of sound. In this analysis, the 
lower order terms are source terms and are assumed 
known (using the procedure provided later). 

Differentiating once more in time provides a new set 
of characteristics with the same velocity as before: 


will be modified such that the boundary condition, 
Eq. (5), is satisfied while simultaneously insuring the 
tangential derivatives of the pressure: 


d a p d a p 
d a r ’ d a ~ 1 Tdr) ’ 


d a p 

dTd a ~ 1 r] 


V(q) : (a <0) 


(22) 


do not change on an inviscid wall. 

We can represent the new pressure, ( )n, as the sum 
of the old, ( )o, and correction, ()c, pressures: 



where 


d(vi) = d ( u(u x (p x - ft-) - u(p xx - E §r))) 
d(v 2 ) -d((u + c)(u x (u x + ^) + {u + c) [u xx + ^)) 

d(v 3 ) =di(u- c)(u x (u x -*£) + («- c)(u xx ~ *£■)) 

( 20 )' 

In Eq. (18) the entropy characteristic, d(v 1 ) does not 
depend on velocity, but the new differentiated charac- 
teristic does depend upon velocity in Eq. (19). And, 
in Eq. (20) the new entropy characteristic depends on 
all flow variable first derivatives. By only correcting 
the pressure derivatives, the Euler system is under- 
specified because the velocity and density derivatives 
are incorrect at the surface. Ideally, all the flow vari- 
able derivatives would be corrected so that d(v 1 ) and 
the outgoing acoustic characteristic are not changed. 
This issue is beyond the scope of this paper since it is 
currently standard practice to only correct the pres- 
sure terms and since physically correct solutions have 
been achieved with it (probably because the required 
corrections in the other terms were small compared to 
the pressure). This will become more important as 
higher order time accuracy is required and the proce- 
dure shown next could be extended for such cases. 

IV. Surface Pressure Correction 
Algorithm 

The basic idea 7 is to modify the normal pressure 
derivatives on a boundary to satisfy the governing 
equations and boundary conditions. One-sided deriva- 
tives at a surface will generally be incorrect since they 
do not account for the walls presence, thus they need 
to be corrected. And, contrary to the assertion of only 
one physical boundary condition being available for 
high order finite-difference schemes, we make use of 
the infinite additional conditions 10 in Eq. (5). 

For a wall with normal vector, if = ( T] x ,T]y ), tangen- 
tial vector, t = (r x ,T y ), and a numerical method with 
order O accuracy, the pressure derivatives: 

^V(a,6 ):(l<a + 6<0) < 21 > 


s, 

a— i,i,0 


(-n.es, o.a - n,CS, o, a + (-1)' M* 


(24) 


3 = 0 


b 


a(a — l) 

2 / a(a- 1 ) 

Mi= J2 (-i) a+I ' 2 

6=0 

A j = 


-.a 2 — 26— i„2b+i 


n 2 y b+l (25) 


-' 0 , 0,0 


—n y (V q ,i,j-i iUiW Co i0i o + Tc'ij- 2,ti,uQ,o,o A ( “'c £ 0 


Ta,pXJ,9 — 
§+! 


/ 1 \ a — 1 V "' IT ( T ' C '?. 0 . 1 )) <,+ ’ T ' ) 

(~ 1 ) 2_j / , Hm,i,a,0X ( c e yi+m+fi) 


(^ 0 , 0 ,o) 


m = 0 i — 0 


(^ 0 , 0 , 0 ) 


(a-(8-2i-C-l) 


(c— 2 m) 


H , 




(a - l)!(m + i)\ 


(a - fd - 2i - C - 1)!(C ~ 2 m)!( 2 (i + m) + p)\i\m\ 


u f 

The updated term, C o,o,a> * n Eq. (24), refers to the 
value of Cq 0 after the lower order pressure deriva- 
tives have been corrected. For example, we first correct 
the 1 st (a = 1 ) order pressure derivatives, (p x ,p y ), us- 
ing Eq. (24). Next, the higher order time derivatives, 

Co 0 2 an d Co, 0,2 are calculated, but using the new val- 
ues of p x and p y . These ’’updated” time derivatives are 
then used to calculate new 2 nd (a = 2 ) order pressure 
derivatives , Pxx,Pxy,Pyy, with Eq. (24). 

This process is repeated until the desired accuracy 
is achieved. A complete derivation of Eq. (24) may be 
found in the appendix. 
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V. Results 

The adequacy of correcting only the pressure deriva- 
tives to achieve high time accuracy was tested by 
scattering the initial Gaussian pulse: 


p{x,y, 0) = 0.001 exp ( — 


ln2[(z- 1 ) 2 + y 2 


.04 


(28) 


with four grid points per unit interval from a flat plate 
at progressively higher order accuracy. The initial uni- 
form flow conditions were p = 1.4, p = 1.0, u = 0, and 
v = 0 

In Fig. 2, pressure contours of the reflected pulse 
are shown for the 3 rd order conditions. The solution is 
stable, but the reflected pressure is slightly amplified. 
The pressure is over-corrected (amplified) because the 
effects of the other uncorrected variables must be ac- 
counted for to maintain the no-flow boundary condi- 
tions, Eq. (5). As the accuracy increases, the number 
and effect of uncorrected variables increases resulting 
in an over-correction of the pressure that leads to nu- 
merical instability, unless some damping mechanism is 
added. 

In Fig. 3, the reflected pulse using the 3 rd order 
conditions are compared with the current practice of 
correcting only the first derivatives, p x and p y using 
3 rd and 1 st order design time accuracy at the surface 
(though both achieve only 1 st order time accuracy as 
mentioned) . The 3 rd order time accuracy cases in Figs. 
3a and 3b, showed little difference here suggesting the 
higher order corrections are small, but the 1 st order 
time advance in Fig. 3c generated spurious waves and 
was unstable. 

The 5 th order example in Fig. 4 shows more clearly 
the importance of correcting the higher derivatives. 
Despite only correcting the pressure, the reflected 
pulse is qualitatively superior to correcting only the 
first order pressure derivatives. The simulation be- 
comes unstable however on the surface after the pulse 
has reflected from it at time t = 0.8. Higher order con- 
ditions result in larger over-corrections of the pressure, 
as expected. Note that the MESA scheme explicitly 
uses the 2 nd order spatial derivatives and the effects of 
high order boundary corrections are directly observ- 
able. Other methods may perform differently for the 
cases shown in Fig. 4b and 4c., though their results 
remain nonphysical. 

As a test of the procedures on curved surfaces with a 
Cartesian grid, the same Gaussian pulse was scattered 
from a sphere of radius, r = .25, with a grid density 
of 128 points per unit interval. In Fig. 5, a 3 rd or- 
der solution is shown. The curvature results in greater 
changes in the velocity and the pressure corrections 
must compensate for that, resulting in a dissipated, 
but clearly reflected pulse and was numerically stable. 
The higher order conditons (> 4) were unstable with- 
out adding artificial dissipation for the circle case. 


As a more realistic example, observe the entropy 
production shown for the 2D stator blade passage in 
Fig. 6 in which the height of the grid indicates mag- 
nitude of the entropy. Only the I s * order pressure 
derivatives are corrected and for this inviscid problem 
the solution should be flat, indicating no change in 
entropy. Here, the curvlinear form of the nonlinear 
Euler equations are solved using a 6 th order compact 
scheme 18 for spatial derivatives and a 4 tft order 5-6 
Low Dispersion and Dissipation Runge-Kutta scheme 
for time marching. 19 Notice the entropy production at 
the suction surface and convected downstream in the 
wake. 

Overall, the advantages of correcting the high or- 
der spatial derivatives is clear, and as shown in Fig. 1, 
the efficiency of high accuracy is evident for demand- 
ing long-term low error tolerance simulations. But the 
current practice of correcting only the pressure is not 
physically correct despite the flow variables satisfying 
the governing equations and boundary conditions at 
the surface. The practice is still nonphysical (e.g., 
produces entropy and vorticity) because the character- 
istics are under-specified and therefore the other flow 
variables must be corrected as well. 

VI. Conclusion 

Most surface boundary methods in use today 
achieve only 1 st order accuracy in time since only the 
no flow condition is imposed. Higher time accuracy 
requires enforcing the highier order time derivatives of 
the no flow condition. In the special case of stationary 
surfaces, the even ordered time derivatives of the no 
flow condition will be automatically satisfied if they 
were correctly specified at the start of the simulation. 

The primary advantage of high time accuracy is 
better resolution and larger time steps at surfaces. 
Achieving high time accuracy requires correcting 
cross-derivatives of the flow at a surface. In principle, 
any numerical method could be used with the bound- 
ary procedure presented here. But a Hermitian MESA 
scheme was used because of its extraordinary resolu- 
tion and potential for efficiency at a surface with high 
time accuracy. The exact performance advantages of 
high time accuracy will depend upon the physical ac- 
curacy required. 

The design time accuracy was increased here by 
correcting only the pressure as is the currently ac- 
cepted l 5t order practice. Despite the apparently high 
complexity in deriving very high order pressure correc- 
tions, a simple process for doing this was found that 
requires only a few lines of code to implement. How- 
ever, a characteristic analysis and numerical testing 
has revealed that it is important to correct the veloc- 
ity and density spatial derivatives as well. Therefore, 
an important next step is to correct all flow derivatives 
to avoid underspecifying the hyperbolic Euler system 
and to extend this analysis to curvilinear coordinates 
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a) n=190 


b) n=690 c) n=1190 

Fig. 5 Third Order Pulse Scattered From Circle 


d) n=1410 


to efficiently resolve curvature effects. 


find the following relationship: 


VII. Appendix: Derivation of Pressure 
Correction 

The pressure correction procedure shown in sec- 
tion IV was derived by assuming only the pressure’s 
spatial derivatives required adjustment to account for 
the presence of the wall. After the corrections are ap- 
plied, the flow variables and their spatial derivatives 
must satisfy the nonlinear Euler equations and the 
boundary condition at the surface. In addition, the 
flow’s spatial derivatives must satisfy the time deriva- 
tives of these equations to achieve high accuracy in 
time. 

Starting with the following boundary condition, 
from Eq. (5), for any a: 


d a (V.n) _ d a u d a v _ nQs 
dt a “ Ux dt a + Uy dt<* ~ 


(29) 


c 

/nrli _ 

'-'0,0, a — 


C u /nU i C'V s~iu 

0,0,0^ 1,0,0: — 1 u 0,0,0'- / 0,l,a- 


+ 


7^7 

^ 0 , 0,0 


r p 

W.O, a -1 


n v — 

'-'0,0, a “ 

~ 0 ^ 1 , 0 , 0-1 + Co,o,oCo,l,a-l + cl 00 C 0 , 1 , a -1 

(32) 

which shows how changes in pressure and velocity 
mixed space and time derivatives will affect the bound- 
ary condition equation. 

The critical and most difficult next step in this 
derivation is expressing the right hand side of Eq. (31) 
in terms of the corrected pressure’s spatial derivatives. 
T his was done by recursively differentiating all six gov- 
erning equation forms with computer algebra until a 
pattern could be discerned. The following pattern was 
observed and inductively verified: 


we want to correct the pressure derivatives so that the 
following is satisfied: 


o 


Co, 0 ,a + C'o.o.a ) + n y ( Co, 0 ,a + Co, 0 ,a J — ^ 0 , 0 , a 

(30) 


c c 

The terms, C o,o,a and C o,o,a> are tiie corrections to 
the time derivatives of the velocities resulting from the 
corrections to the spatial derivatives of the pressure. 
We need to determine how changes in the pressure 
derivatives will effect the velocity time derivatives so 
that Eq. (30) is satisfied up to some specified order a. 

This relationship can be found using the governing 
equations, Eqs. (1, 2, 3, 4) and the following alternate 
form of the momentum equations: 16 


Cl, 0, a-1 

= C P a -jjfiT a ,i,i,u,v 

3= 0 

C 

/'-m 

^0,1, a-1 

^ Q 

3=0 

CUa-l 

“ c 

= T. Ca-j'j'oTa,l,j-l,u,v 

3=0 

Cl i, a _i 

™ Q 

3=0 

r p 

W,0,a-1 

“ C 

= C'a-j,j, 0 T a, 0 ,j,u,v 

3=0 

c p 

0,l,a— 1 

“ c 

3=0 


Su _ _ ( U §u , du , 1|£ 

dt l ^ dx dy p dx 

§v _ _ L,dv , dv , 1 9£ 

dt ~ \ U dx ^ V dy ^ p dy 


(31) 


where is defined in Eq. (27). 

For each value of a, we can now correct the following 
pressure derivatives: 


From this alternate momentum equation form we 


Ca-nnoVi=0,l, 


(34) 
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Fig. 6 Entropy Produced on Surface of Stator Blade On Curvilinear Grid 
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by subsituting Eq. (33) in to Eq. (32) and solving 
boundary condition Eq. (30): 


Ca-j'jfiAj — CQ' Q a n x Co, 0 : a n yC 0 , 0 , a ( 35 ) 


j = 0 


with Aj defined in Eq. (26). 

In addition, since we do not want corrections of the 
pressure derivative to effect the tangential derivatives, 
we impose, for each a the following a conditions : 


d a P \ 
dr 01 ~i dpi ) 


correction 


= 0 Vj:i = 0 ,l,.. 


a — 1 
(36) 


which can be written as a condition of only the Carte- 
sian x and y pressure derivatives as: 


where matrix F equals matrix E with column number 
i replace by vector D. 

The determinant of a matrix may be solved using 
cofactor expansion along a column. 1 ' So we can take 
advantage of all the zeros in F by expanding along 
column i to efficiently evaluate the determinant of the 
matrix: 


co factor 

det(F) — (Cqq 0 — Tl x Co,0,a — n yC 0,0, a) 

(40) 

where Mi is the minor of the element in the first row 
and i th column of matrix F. 

Similarly we can evaluate the determinant of matrix 
E by cofactor expansion along the top row: 


d a p 


dT a ~ j dpi / correction. 

= ( T *£;+Tyi l) J (vx£+Py-^y P 

= Ca-a,a,0-E?a,j 

a=0 

= o 


(37) 


with 



6=0 


a- j 
a — b 



(_l)(i>-a)( rlx )0-26+a)( n ^^(a-i+26-a) 


since t x = p y and r y = -p x . 

For each value of a, we can write the a+1 conditions 
in Eqs. (35,37) in the following matrix form: 






c„ 1 

Ao 

A\ ... 

A a 


a, 0,0 

#o,o 

#1,0 

#a,0 


% 

^ a— 1,1,0 
C 

#o,i 

#1,1 . . . 

#a, 1 


#0,2 

#1,2 

#a, 2 


^a- 2,2,0 

#0,a-l 

Bl t a— 1 



c p 


n xCo, 0 ,a n yCo ,0 ,q 
0 
0 
0 


u u 


(38) 


D 


Since the right hand side of this matrix system is 
zero, except the first term, we can efficiently solve this 
system with Cramer’s rule for the i th element of vector 
C: 

det (F) 

* det(E) 


det(E) = J2A j (-l) j M j 

j= o 

(-i)°K+"^ ( Yl n *- in v A 

i = 0 


(41) 


where Aj is defined in Eq. (26). The expression for 
Mj, from Eq. (25), was derived by observation with 
the help of computer algebra and verified inductively. 

With Mj known, we can quickly apply Cramer’s rule 
to solve Eq. (39) for the unknown pressure correction 
vector C in Eq. (38). This simple procedure is repeated 
for a — 1, 2, . . . , O until all higher order pressure spa- 
tial derivatives are corrected. 
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